package com.asa.numanaly;

import com.asa.HanShu;
import com.asa.HanShuN;
import com.asa.hutils.MathHUtils;
import com.asa.utils.MathUtils;

/**
 * 非线性方程组的数值解
 * 这个求解必须是要给定初始值的，而且是可能存在多个零点，而这里只能找初始值附近的解
 * 
 * @author Administrator
 *
 */
public class J {

	/**
	 * 方程组的Newton方法
	 * 
	 * @return
	 */
	public static double shuzhijie() {

		return 0;
	}

	public static void main(String[] args) throws Exception {

		HanShuN hs1 = new HanShuN() {

			@Override
			public double hanshu(double t, double[] x) {
				// TODO Auto-generated method stub
				return 0;
			}

			@Override
			public double hanshu(double[] x) {
				// TODO Auto-generated method stub
				return 3 * x[0] - Math.cos(x[1] * x[2]) - 0.5;
			}
		};
		HanShuN hs2 = new HanShuN() {

			@Override
			public double hanshu(double t, double[] x) {
				// TODO Auto-generated method stub
				return 0;
			}

			@Override
			public double hanshu(double[] x) {
				// TODO Auto-generated method stub
				return x[0]*x[0] - 81.0*(x[1]+0.1)*(x[1]+0.1) +Math.sin(x[2])+1.06;
			}
		};
		HanShuN hs3 = new HanShuN() {

			@Override
			public double hanshu(double t, double[] x) {
				// TODO Auto-generated method stub
				return 0;
			}

			@Override
			public double hanshu(double[] x) {
				// TODO Auto-generated method stub
				return Math.exp(-x[0]*x[1])+20*x[2]+(10*Math.PI-3)/3;
			}
		};
		double[] x = {0.1,0.1,-0.1};
		HanShuN[] basehs = new HanShuN[3];
		
		basehs[0] = hs1;
		basehs[1] = hs2;
		basehs[2] = hs3;
		int k = 5;
//		fangchengzu2(x, basehs,k,-1);//不想提前出来，就把精度设为负数即可
		
		
		fangchengzu1(x, basehs, k, -1);

	}

	/**
	 * 牛顿法
	 * 我感觉这个比拟牛顿法简单，但我是先完成的拟牛顿法，因为那个算法给了例子，我根据例子反推的算法
	 * 所以很多东西就是通用的
	 * 
	 * @param x
	 * @param basehs
	 * @param k
	 * @param TOL
	 * @return
	 * @throws Exception
	 */
	private static double[] fangchengzu1(double[] x, HanShuN[] basehs,int k,double TOL) throws Exception {
		
		double[] next_x = new double[x.length];
		for (int i = 0; i < x.length; i++) {
			next_x[i] = x[i];
		}
		
		
		
	
		
		
		for (int i = 0; i < k; i++) {
			double[][] jakebi = jakebi(basehs, next_x);
			double[] yuanfunc_value = yuanfunc_value(basehs, next_x);
			double[] s = F.fangchengzu1(jakebi, yuanfunc_value);

			
			for (int j = 0; j < s.length; j++) {
				next_x[j] = next_x[j] - s[j];
				System.out.println(next_x[j]);
			}
			double TOL_ = 0;
			for (int j = 0; j < s.length; j++) {
				TOL_ += s[j]*s[j];
			}
			TOL_ = Math.pow(TOL_, 0.5);;
			System.out.println("精度是："+TOL_);
			if (TOL_<TOL) {
				return next_x;
			}
			
		}
		
		
		
		return next_x;
		
	}

		
	
	
	
	
	/**
	 * 拟牛顿法求非线性方程组的解
	 * 
	 * p555??	Broyden，这个是算法的名字，但提供出来的算法太简洁了，导致阅读有些困难
	 * 所以我实际上是根据P556给出的例子，自己完成的算法步骤
	 * 这里的这个是三维的，我不确定这个实现是否已经支持不同于三维的求解，用到时，需要确认
	 * @param x	初始的x值，这个其实和求方程的零点一样的，用牛顿法来求，就是要初始值，然后用导数来靠近零点
	 * @param basehs	函数方程组
	 * @param k	最大迭代次数
	 * @param TOL 提前跳出循环的精度		如果不想提前出来，就把精度设为负数即可
	 * @throws Exception
	 */
	private static double[] fangchengzu2(double[] x, HanShuN[] basehs,int k,double TOL) throws Exception {
		double[][] move2nextAn = move2nextA0(basehs, x);
		
		double[] yuanfunc_value = yuanfunc_value(basehs, x);

		
		double[][] next_xands = move2nextX(x, yuanfunc_value, move2nextAn);
		double[] next_x = next_xands[0];
		double[] s = next_xands[0];

		
		for (int i = 1; i <= k; i++) {
			move2nextAn = move2nextAn(basehs, yuanfunc_value, s, next_x,move2nextAn);
			yuanfunc_value = yuanfunc_value(basehs, next_x);
			
			next_xands = move2nextX(next_x, yuanfunc_value, move2nextAn);
			next_x = next_xands[0];
			s = next_xands[1];
			System.out.println("每次的x的值+"+i);
			for (int j = 0; j < next_x.length; j++) {
				System.out.println(next_x[j]);
			}
			
			double TOL_ = 0;
			for (int j = 0; j < s.length; j++) {
				TOL_ += s[j]*s[j];
			}
			TOL_ = Math.pow(TOL_, 0.5);;
			System.out.println("精度是："+TOL_);
			if (TOL_<TOL) {
				return next_x;
			}
		
			
		}
		return next_x;
	}
	
	
	private static double[][] move2nextX(double[] xpre,double[] yuanfunc_value,double[][] A ){
		double[][] fx = new double[1][];
		fx[0] = yuanfunc_value;
		double[][] chenfa = MathHUtils.chenfa(fx, A);

		double[] next_x = new double[xpre.length];
		double[] s = new double[xpre.length];

		for (int i = 0; i < next_x.length; i++) {
			next_x[i] = xpre[i]-chenfa[0][i];
			s[i] = -chenfa[0][i];
			System.out.println(s[i]+"     s...    ");
			
			

		}
		
	
		double[][] result = new double[2][];
		result[0] = next_x;
		result[1] = s;
		return result;
		
	}

	
	
	
	
	private static double[][] move2nextAn(final HanShuN[] basehs,double[] fx,double[] s, final double[] next_x,double[][] A ){
		
		double[] yuanfunc_value = yuanfunc_value(basehs, next_x); 
		
		double[] y1 = new double[yuanfunc_value.length];
		
		for (int i = 0; i < y1.length; i++) {
			y1[i] = yuanfunc_value[i] - fx[i];
			System.out.println(y1[i]+"=========y1");
		}
		double suanzi = 0;
		double[] suanzizu = new double[s.length];
		
		for (int i = 0; i < suanzizu.length; i++) {
			
			for (int j = 0; j < suanzizu.length; j++) {
				suanzizu[i] += s[j]*A[i][j];
			}
		}
		for (int i = 0; i < suanzizu.length; i++) {
			suanzi+=suanzizu[i]*y1[i];
		}
		System.out.println("suanzi      "+suanzi);
		
		//拿到算子后，开始计算下一个A逆
		double[] A_y1 = new double[s.length];
		
		for (int i = 0; i < A_y1.length; i++) {
			for (int j = 0; j < suanzizu.length; j++) {
				A_y1[i] += A[i][j]*y1[j];
			}
		}
		double[] s__ = new double[s.length];
		
		for (int i = 0; i < s__.length; i++) {
			s__[i] = s[i] - A_y1[i];
		}
		double[][] A_pre = new double[s.length][s.length];
		for (int i = 0; i < A_pre.length; i++) {
			for (int j = 0; j < A_pre.length; j++) {
				A_pre[i][j] += s__[i]* suanzizu[j];
			}
		}
		double[][] chenfa = MathHUtils.chenfa(A_pre, 1.0/suanzi);
		double[][] jiafa = MathHUtils.jiafa(A, chenfa);
		for (int i = 0; i < jiafa.length; i++) {
			for (int j = 0; j < jiafa.length; j++) {
				System.out.print(jiafa[i][j]+"     ");
			}
			System.out.println();
		}
		
		
		return jiafa;

	}

	
	
	
	
	
	
	
	
	

	private static double[][] move2nextA0(final HanShuN[] basehs, final double[] x)
			throws Exception {
//		double[][] Jx = new double[3][3];
		
		
		double[][] jkbx = jakebi(basehs, x);
		
//		double[] fx = new double[3];//原函数的值
//		yuanfunc_value(basehs, x, fx);
//		for (int i = 0; i < fx.length; i++) {
//			System.out.println(fx[i]);
//		}
		
		
		
		
		
		for (int i = 0; i < jkbx.length; i++) {//这个就是一次的雅克比矩阵
			for (int j = 0; j < jkbx.length; j++) {
				System.out.print(jkbx[i][j]+"      "  );
			}
			System.out.println();
		}
		
		System.out.println("=============");
		for (int i = 0; i < x.length; i++) {
			System.out.println(x[i]);
		}
		
		System.out.println("=============");

//		MathUtils.dao(x0, hanShu)
		
		
		double[][] fangchengzu000 = F.fangchengzu000(jkbx);
//		double[][] zhuanzhi = MathHUtils.zhuanzhi(fangchengzu000);
//		
		System.out.println("第一次求解出来的逆：");//这个行列是不需要转置的，但和书上比起来是转置了一下，但不需要
		for (int i = 0; i < fangchengzu000.length; i++) {//雅克比求逆
			for (int j = 0; j < fangchengzu000[i].length; j++) {
				System.out.print(fangchengzu000[i][j]+"      "  );
			}
			System.out.println();
		}
		
		
		
		
		return fangchengzu000;
		
		
		
	}
	
	
	/**
	 * 雅克比矩阵的值
	 * @param basehs
	 * @param x
	 * @return
	 */
	private static double[][] jakebi(final HanShuN[] basehs, final double[] x) {
		HanShu[][] jkb = new HanShu[3][3];

		
		for (int i = 0; i < basehs.length; i++) {
			for (int j = 0; j < basehs.length; j++) {
				final int i2 = i;
				final int j2 = j;
				jkb[i][j] = new HanShu() {
					
					@Override
					public double hanshu(double t, double x) {
						// TODO Auto-generated method stub
						return 0;
					}
					
					@Override
					public double hanshu(double x0) {
						// TODO Auto-generated method stub
						double[] copyx = new double[x.length];
						for (int i = 0; i < copyx.length; i++) {
							copyx[i] = x[i];
						}
						copyx[j2]=x0;
						return basehs[i2].hanshu(copyx);
					}
				};
				
				
			}
		}
		
		
		double[][] jkbx = new double[3][3];

		for (int i = 0; i < basehs.length; i++) {
			for (int j = 0; j < basehs.length; j++) {
				
				jkbx[i][j] = MathUtils.dao(x[j], jkb[i][j]);
				
			}
		}
		return jkbx;
	}

	/**
	 * 原方程的带入x的值
	 * @param basehs
	 * @param x
	 * @return
	 */
	private static double[] yuanfunc_value(final HanShuN[] basehs, final double[] x) {
		double[] fx = new double[basehs.length];
		for (int i = 0; i < basehs.length; i++) {
			fx[i] = basehs[i].hanshu(x);
		}
		return fx;
	}



}
